Destruction of Anderson localization by a weak nonlinearity 
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We study numerically a spreading of an initially localized wave packet in a one-dimensional discrete 
nonlinear Schrodinger lattice with disorder. We demonstrate that above a certain critical strength 
of nonlinearity the Anderson localization is destroyed and an unlimited subdiffusive spreading of 
the field along the lattice occurs. The second moment grows with time oc t a , with the exponent a 
being in the range 0.3 — 0.4. For small nonlinearities the distribution remains localized in a way 
similar to the linear case. 
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The Anderson localization [l[ has been originally dis- 
cussed for electron propagation in a disordered potential. 
Nowadays, an impressive technological progress in exper- 
iments with cold atoms allows one to create a disordered 
quasi-lD potential by laser fields and to detect signatures 
of localization of a Bose-Einstein condensate (BEC) in 
presence of disorder 0, H, 0, 0, @ ■ An interesting new as- 
pect in such systems is an importance of nonlinear effects 
since in a good approximation the evolution of BEC can 
be described by the nonlinear Gross-Pitacvskii equation 
(see e.g. 0). An interplay of disorder, localization, and 
nonlinearity appears also in other physical systems like 
wave propagation in nonlinear disordered media (see e.g. 
d, Q), chains of nonlinear oscillators (see e.g. [Ic]]) with 
randomly distributed frequencies, and models of quan- 
tum chaos with a kicked soliton and a kicked BEC 

[II El El. 

We focus here on the discrete Anderson nonlinear 
Schrodinger equation (DANSE) 



ih 



dt 



E n 1p„ + (3\ V>„ | Ipn + V(l/J n+1 + V„-l) , (1) 



where (3 characterizes nonlinearity, V is hopping matrix 
element, on-site energies are randomly distributed in the 
range — W/2 < E„ < W/2, and the total probability is 
normalized to unity \ tp n \ =1- For (3 = and weak 
disorder all eigenstates are exponentially localized with 
the localization length I rs 96(V/W) 2 at the center of 
the energy band [l9(. Hereafter we set for convenience 
H = V = 1, thus the energy coincides with the frequency. 

For nonlinear equation ([T]) we consider the following 
problem: how an initially localized field | yj n (0) | 2 = S n fl 
is spreading? In the linear case the spreading saturates 
after excitation of all linear modes that have significant 
values at n = 0. The same process of "initial excita- 
tion" of modes happens in the nonlinear case as well, 
this initial stage of spreading has been analyzed recently 



in refs. [20j, |2l| and is now well understood. However, 
a behavior at large time scales remains less clear. The 
results presented in support the view of eventual ex- 
ponential localization of the field. We demonstrate below 



that the spreading is unlimited, however it is rather slow 
- subdiffusive. 

The basic idea is to use the equivalence between the 
Anderson localization and the localization ofquasienergy 
eigenstates in a kicked quantum rotator 0, El- In 
the latter model the case of quantum chaos with nonlin- 
earity has been considered analytically and numerically 
in [13, EH and it has been shown that above a certain 
nonlinearity level, nonlinear phase shifts lead to a com- 
plete delocalization with a subdiffusive spreading over all 
states [l7| ■ Furthermore it has been argued that the same 
situation should appear for the DANSE ([I]). 

We first apply the theoretical arguments of paper [l7| 
to model ([1]) , and then perform large scale numerical sim- 
ulations of a wave packet spreading on a time scale which 
is by 5-6 orders of magnitude larger compared to that in 
jifit . Our results are in general consistent with the theory 
developed for the quantum chaos model [17[ that predicts 
unlimited subdiffusive spreading. 

At first glance, the effect of nonlinearity seems to be 
vanishing in the limit of a broadly spread distribution. 
Indeed, if the field is spread over An sites, then due to 
the conservation of the total probability in Eq. ([T]) the 
field is small | ip n | 2 ~ 1/An and correspondingly small 
arc the nonlinear effects. However, one should compare 
the nonlinear frequency shift Su> ~ (3 \ ip n | 2 with the 
characteristic distance Aw between frequencies of excited 
modes (the latter are the exponentially localized modes 
of the linear disordered lattice). As the number of these 
modes is proportional to An, the distance between the 
frequencies obeys Aw ~ 1/An, and the relative nonlin- 
ear frequency shift Slo/Auj ~ (3 is independent of the 
width of the field distribution but is proportional to the 
nonlinearity parameter (3. This means that the effect of 
nonlinearity does not qualitatively depend on the width 
of field distribution: if for some (3 > (3 C the field is 
chaotic, chaos remains while spreading, and no transi- 
tion to regularity that blocks spreading is expected; for 
small nonlinearities (3 < (3 C the dynamics which is nearly 
regular (KAM regime) for localized distributions remains 
as such for all times. However, quantitative difference 



2 



can occur and the spreading can slow down for wide dis- 
tributions. Again, to roughly estimate this effect, we 



can adopt here the arguments of [17[. In the basis of 
linear localized modes, the evolution of the amplitudes 
C m of these modes is due to their nonlinear coupling, 
i.e. C m ~ f3C mi C m2 C m3 . Assuming randomness of the 
phases, we can estimate the rate of excitation of a newly 
involved mode as ~| C l/(An) 3 . On the other hand, 
excitation of a new mode is none other than diffusive 
spreading of the field, thus |(An) 2 ~ l/(An) 3 . Solu- 



tion of this equation yields subdiffusivc spreading 
(An) 2 cx t 2 '^ . 



(2) 



For numerical simulations we used the operator split- 
ting integration scheme for the time evolution given by 
©: ip n (t + At) = R§if>„(t), where R = exp(-i(E n + 
(3\ ipn | 2 ) At) is local and S = exp(— 2 At cos#) is nontriv- 
ial because 9 is the operator conjugated to n = —id/d6. 
This kick-like integration scheme is unitary and preserves 
the total probability. In addition it introduces high 
harmonics with frequencies ui m = mlit j At and integer 
to. However, at small At these frequencies arc signif- 
icantly larger than the system energy band width B: 
lo\ = 2n/At 3> B = (4 + W) and their average effect 
should be exponentially small. We have chosen At = 0.1 
that gives high frequency oscillations of total energy on 
a percent level. A further decrease of At by an order 
of magnitude does not affect the average behavior of the 
field spreading. 

We used two discrete implementations for the evolu- 
tion operator S of the linear Schrodinger equation. In 
the first one we represent S as a band matrix whose ele- 
ments are Bessel functions. At At = 0.1 we keep Besscl 
functions J m (2Ai) with | to |< 10 that preserves prob- 
ability on one integration step with the accuracy better 
than 10 -16 , whereas after time t = 10 8 the probability is 
preserved with accuracy better than 10~ 7 . For one disor- 
der realization such a run with the total number of sites 
N = 2001 and t = 10 s takes about 6000 mins of CPU 
on a iGHz- workstation. In another implementation we 
used the unitary Crank-Nicholson scheme (22|. The re- 
sults obtained by both methods are very similar; below 
mainly those from the first method are presented because 
of its better efficiency. 

To characterize the wave packet spreading over the lat- 
tice sites we compute its average squared width, i.e. the 
second moment ((An) 2 ) = a(t) = J2n( n ~ ( n )) 2 \ V^iW | 2 - 
The averaging over disorder realizations was performed 
for the logarithm of this quantity, i.e. for logo - . The de- 
pendence of the averaged a on time t is shown in Fig. [1] 
for a moderate nonlinearity f3 — 1 and disorder strengths 
W = 2 and W = 4. It clearly shows a subdiffusivc 
spreading 




a(t) oc t a 



(3) 
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FIG. 1: (color online) The dependence of the second moment 
a of the probability distribution w n on time t for two disorder 
strengths W — 2 (top full red/gray curve) and W = 4 (bottom 
full black/blue curve) for /3 = 1; dotted curves of same color 
show data at /3 = (for large t only the average value is shown 
by a horizontal line) . The values of log 10 a are averaged over 8 
disorder realizations and over time intervals A(log 10 t) ~ 0.1. 
Dashed lines show the numerical fits log 10 a = alog 10 t + 77 
with a = 0.344±0.003, n = 1.76±0.02 (done for 3 < log 10 t < 
8 for W = 2) and a = 0.306 ±0.002, n = 0.94 ±0.01 (done for 
2 < log 10 t < 8 for W = 4). The full straight line shows the 
slope 01 = 0.4. Here and below the logarithms are decimal. 



which continues without any sign of saturation up to ex- 
tremely large time t max = 10 8 . At t max the variance 
a becomes by two orders of magnitude larger than its 
saturation value at [3 — 0. Yet the initial spreading at 
t < I for j3 = 1 is similar to the linear case (3 — in 
agreement with [2l| . The exponent a was determined by 
a fit in the time interval t$ < t < t ma x, where to is the 
characteristic time at which the linear spreading ends: 
a((3 = 0;t) < a(f3, to). The fits are shown in Fig. Q] and 
the fit values are given in the caption to Fig. [TJ The sta- 
tistical error in the value of a is rather small due to a 
large time interval, however a values for individual real- 
izations fluctuate rather strongly varying in the interval 
0.32 < a < 0.39 and 0.28 < a < 0.41 with the standard 
deviation errors Aa/a = 0.026 and 0.045, for IT = 2 and 
4 respectively. There are also certain time variations, e.g. 
for 10 5 < t < 10 8 the fits give a = 0.375 and 0.319 for 
W = 2 and 4 respectively. In spite of these variations the 
value of a differs visibly from the theoretically expected 
value a = 0.4 (Ref. [lj and Eq. j2|)). Wc will return to 
the discussion of this deviation later. 

A more detailed characterization of the field spreading 
is given by the probability distribution w n = | ip n (t) | 2 . 
We show log 10 w n averaged over disorder realizations in 
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FIG. 2: (color online) Probability distribution w n over lattice 
sites n at W = 4 for (3 = 1, t = 10 8 (top blue/solid curve) 
and t = 10 (middle red/gray curve); /3 = 0, £ = 10 5 (bottom 
black curve; the order of the curves is given at n = 500). 
At (3 = a fit lnui„ = — (-y|n| + x) gives 7 w 0.3, x ~ 4. 
The values of log 10 w n are averaged over the same disorder 
realizations as in Fig. [T] 



Figs. EH at t = 10 5 and 10 s . For W = 4 the tails of the 
probability distribution drop down to enormously small 
values ~ 10 -130 that can be reached due to our inte- 
gration scheme which works efficiently up to very small 
absolute values of probability. The tails of the distribu- 
tion w n drop exponentially with the same slope as for 
the linear case (3 = which is also shown; the decay rate 
7 ps 0.30 is close to the theoretical value 2/1 ss 0.33. An- 
other notable feature of w n is a flat distribution, chapeau, 
centered near the initially populated site n = 0. Inside 
the chapeau the sites are populated in an approximately 
homogeneous way, and hence its width is essentially de- 
termined by the second moment u(t). For W = 2 the 
decay rate 7 for /3 = drops approximately by a factor 
ps 4 compared to the case W = 4, in agreement with 
the theoretical expression for /. Due to a larger value of 
I, the spreading over the lattice sites is broader and a 
non-exponential shape of the distribution w n at t = 10 8 
is more visible. At shorter times t = 10 5 the tails of the 
distribution are very similar to those in the linear case 
[3 = 0. 

The dependence of a(t) on nonlincarity (3 is shown in 
Fig. U] for W = 4. One can clearly see that for (3 = 0.03 
the growth of a is stopped completely, the probability 
distribution w„ is close to exponential localization like 
for (3 = (see Fig. 5). For (3 = 0.1 there is still a very 
slow increase of a with time, which is however hardly 
distinguishable from a saturation. This value of (3 is pre- 
sumably close to a critical one, at which the unlimited 
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FIG. 3: (color online) Same as in Fig. [2] but with W = 2. 
At p = a fit lnw n = — (7|n| + x) gives 7 « 0.06, x ~ 
—3. The values of In w n are averaged over the same disorder 
realizations as in Fig. [T] 

spreading sets on. A similar transition occurs for W = 2. 
These data show that a derealization transition takes 
place at a certain critical nonlincarity (3 C « 0.1. We note 
that the qualitative and quantitative features of the dy- 
namics seem to be independent of the sign of (3: the 
spreading for (3 = — 1 is similar to that for (3=1 (we 
have checked this also for (3 = ±2). 
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FIG. 4: (color online) Dependence of the second moment a of 
probability distribution w n on time t for different strengths 
of nonlinearity (3 = 1, 0.1, 0.03, at W = 4 (curves from top 
to bottom at log 10 £ = 4.5, respectively). Data are shown for 
one particular disorder realization. 
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FIG. 5: (color online) Probability distribution w n over lattice 
sites n for the same disorder realization as in Fig. [4] for W = 
4, t = 10 s and for different nonlinearities (3 = 1,0.1,0.03,0 
(from top to bottom; for clarity 3 bottom curves have ad- 
ditional vertical shifts -10, -20 and -30 compared to the top 
one). 



The obtained results show close similarities between 
the DANSE model and the nonlinear kicked rotator stud- 
ied in 



171 ]. In both cases for > /3 C a subdiffusive 



spreading over the lattice continues up to enormously 
long times. In both models the probability distribution 
has a chapeau with approximately homogeneous proba- 
bility distribution. Outside of the chapeau the probabil- 
ity drops exponentially. The width of the chapeau grows 
in a subdiffusive way and the exponent a of this growth 
is approximately the same for both models. The question 
about the exact value of the exponent remains open. It 
is possible that at W = 4 the localization length is rela- 
tively short and there are deviations from the theoretical 
value a = 0.4 given in [l7j]. Longer computations with a 
better statistical averaging are needed to determine the 
exact value of a; the latter may also depend on the pa- 
rameters j3 and W. In particular, slower diffusion might 
be due to inhomogeneities of the effective mode-to-mode 
hopping rates, which arc more pronounced for smaller lo- 
calization lengths, i.e. for larger disorder W. At the same 
time the obtained numerical data clearly show the exis- 
tence of unlimited spreading over the lattice for (3 > (3 C . 
Indeed, for the data of Fig. Q] at W — 2 the nonlinear fre- 
quency shift 6(j w j3w n « ft / \/ <j(t max ) « 0.006 is much 
smaller than the typical level spacing between localized 
states 8v m B jl pa 0.25. The same is true for W = 4. 
This means that in simulations we have reached the time 
scales with apparently asymptotic behavior. 

In conclusion, we have demonstrated that in a one- 
dimensional nonlinear disordered lattice the Anderson lo- 
calization is destroyed and the field spreads subdiffusively 



far beyond the linear localization range. This effect ap- 
pears to have a threshold in the nonlinearity coefficient, 
although the transition might be not perfect, as even for 
small nonlinearities an extremely slow spreading of the 
field due to Arnold diffusion mechanism is not excluded. 
It appears also promising to look for other manifesta- 
tions of nonlinearity-induced destruction of localization, 
e.g. in the scattering problem [23[. 



We thank S.Aubry and S. Flach for useful discussions. 
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